#The following code allows one to replicate all weekly analyses pertaining to Somalia
#it uses the following files: weekly_cyclone_ts.dta and pests.r

#load necessary libraries
library(foreign)
library(MASS)

#source in pests R code
source("pests.r")

#set random number seed
set.seed(10)

#read-in dataset
cyclonedata <- read.dta(file = "weekly_cyclone_ts.dta")

#attach data
attach(cyclonedata)

############
##Select p##
############

# PAR(1) Model
PAR1 <- Parp(som_rebcit_acled ~ 1, p=1)
print(PAR1)
# PAR(2) Model
PAR2 <- Parp(som_rebcit_acled ~ 1, p=2)
print(PAR2)
# PAR(3) Model
PAR3 <- Parp(som_rebcit_acled ~ 1, p=3)
print(PAR3)
# PAR(4) Model
PAR4 <- Parp(som_rebcit_acled ~ 1, p=4)
print(PAR4)
# PAR(5) Model
PAR5 <- Parp(som_rebcit_acled ~ 1, p=5)
print(PAR5)
# PAR(6) Model
PAR6<- Parp(som_rebcit_acled ~ 1, p=6)
print(PAR6)

#########################################
##Select proper exogenous variable lags##
#########################################

#small
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+som_govreb_acled, p=1)
print(Cyclone.PAR1)
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+som_govreb_acled+lagsom_govreb_acled, p=1)
print(Cyclone.PAR1)
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+som_govreb_acled+lagsom_govreb_acled+lag2som_govreb_acled, p=1)
print(Cyclone.PAR1)
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+lagfullcyclone+som_govreb_acled, p=1)
print(Cyclone.PAR1)
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+lagfullcyclone+som_govreb_acled+lagsom_govreb_acled, p=1)
print(Cyclone.PAR1)
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+lagfullcyclone+som_govreb_acled+lagsom_govreb_acled+lag2som_govreb_acled, p=1)
print(Cyclone.PAR1)
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+lagfullcyclone+lag2fullcyclone+som_govreb_acled, p=1)
print(Cyclone.PAR1)
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+lagfullcyclone+lag2fullcyclone+som_govreb_acled+lagsom_govreb_acled, p=1)
print(Cyclone.PAR1)
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+lagfullcyclone+lag2fullcyclone+som_govreb_acled+lagsom_govreb_acled+lag2som_govreb_acled, p=1)
print(Cyclone.PAR1)

#medium
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+som_govreb_acled+lagsom_govreb_acled+som_rebgov_acled, p=1)
print(Cyclone.PAR1)
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+som_govreb_acled+lagsom_govreb_acled+som_rebgov_acled+lagsom_rebgov_acled, p=1)
print(Cyclone.PAR1)
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+som_govreb_acled+lagsom_govreb_acled+som_rebgov_acled+lagsom_rebgov_acled+lag2som_rebgov_acled, p=1)
print(Cyclone.PAR1)

#large
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+som_govreb_acled+lagsom_govreb_acled+som_rebgov_acled+som_govcit_acled, p=1)
print(Cyclone.PAR1)
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+som_govreb_acled+lagsom_govreb_acled+som_rebgov_acled+som_govcit_acled+lagsom_govcit_acled, p=1)
print(Cyclone.PAR1)
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+som_govreb_acled+lagsom_govreb_acled+som_rebgov_acled+som_govcit_acled+lagsom_govcit_acled+lag2som_govcit_acled, p=1)
print(Cyclone.PAR1)

#extra large
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+som_govreb_acled+lagsom_govreb_acled+som_rebgov_acled+som_govcit_acled+lagsom_govcit_acled+som_ncd+russian, p=1)
print(Cyclone.PAR1)
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+som_govreb_acled+lagsom_govreb_acled+som_rebgov_acled+som_govcit_acled+lagsom_govcit_acled+som_ncd+lagsom_ncd+russian, p=1)
print(Cyclone.PAR1)
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+som_govreb_acled+lagsom_govreb_acled+som_rebgov_acled+som_govcit_acled+lagsom_govcit_acled+som_ncd+lagsom_ncd+lag2som_ncd+russian, p=1)
print(Cyclone.PAR1)
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+som_govreb_acled+lagsom_govreb_acled+som_rebgov_acled+som_govcit_acled+lagsom_govcit_acled+som_ncd+russian+lagrussian, p=1)
print(Cyclone.PAR1)
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+som_govreb_acled+lagsom_govreb_acled+som_rebgov_acled+som_govcit_acled+lagsom_govcit_acled+som_ncd+lagsom_ncd+russian+lagrussian, p=1)
print(Cyclone.PAR1)
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+som_govreb_acled+lagsom_govreb_acled+som_rebgov_acled+som_govcit_acled+lagsom_govcit_acled+som_ncd+lagsom_ncd+lag2som_ncd+russian+lagrussian, p=1)
print(Cyclone.PAR1)
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+som_govreb_acled+lagsom_govreb_acled+som_rebgov_acled+som_govcit_acled+lagsom_govcit_acled+som_ncd+russian+lagrussian+lag2russian, p=1)
print(Cyclone.PAR1)
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+som_govreb_acled+lagsom_govreb_acled+som_rebgov_acled+som_govcit_acled+lagsom_govcit_acled+som_ncd+lagsom_ncd+russian+lagrussian+lag2russian, p=1)
print(Cyclone.PAR1)
Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+som_govreb_acled+lagsom_govreb_acled+som_rebgov_acled+som_govcit_acled+lagsom_govcit_acled+som_ncd+lagsom_ncd+lag2som_ncd+russian+lagrussian+lag2russian, p=1)
print(Cyclone.PAR1)


##############
##Table A.11##
##############

Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+som_govreb_acled+lagsom_govreb_acled, p=1)
print(Cyclone.PAR1)

Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+som_govreb_acled+lagsom_govreb_acled+som_rebgov_acled, p=1)
print(Cyclone.PAR1)

Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+som_govreb_acled+lagsom_govreb_acled+som_rebgov_acled+som_govcit_acled+lagsom_govcit_acled, p=1)
print(Cyclone.PAR1)

Cyclone.PAR1 <- Parp(som_rebcit_acled ~ fullcyclone+som_govreb_acled+lagsom_govreb_acled+som_rebgov_acled+som_govcit_acled+lagsom_govcit_acled+som_ncd+russian, p=1)
print(Cyclone.PAR1)


